Parallel-tempering cluster algorithm for computer simulations of critical phenomena 
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In finite-size scaling analyses of Monte Carlo simulations of second-order phase transitions one 
often needs an extended temperature range around the critical point. By combining the parallel 
tempering algorithm with cluster updates and an adaptive routine to find the temperature window 
of interest, we introduce a flexible and powerful method for systematic investigations of critical 
phenomena. As a result, we gain one to two orders of magnitude in the performance for 2D and 
3D Ising models in comparison with the recently proposed Wang-Landau recursion for cluster algo- 
rithms based on the multibondic algorithm, which is already a great improvement over the standard 
multicanonical variant. 



PACS numbers: 05.10.Ln,02.70.Uu,64.60.Cn 



I. INTRODUCTION 

While much attention has been paid in the past to 
simulations of first-order phase transitions and systems 
with rugged free-energy landscapes in generalized en- 
sembles (umbrella, multicanonical, Wang-Landau, paral- 
lel/simulated tempering sampling) [l| , the merits of this 
non-Boltzmann sampling approach also for simulation 
studies of critical phenomena have been pointed out only 
recently. In Ref. [4I , Berg and one of the authors com- 
bined multibondic sampling J3| with the Wang-Landau 
recursion [J| to cover the complete "desired" critical tem- 
perature window in a single simulation for each lattice 
size, where the "desired" range derives from a careful 
finite-size scaling (FSS) analysis of all relevant observ- 
ables. 

Recent developments in the field of Graphic Process- 
ing Units (GPUs) make it possible to have access to a 
massively parallel computing solution at a low cost. To 
use these devices in a most efficient way new parallelized 
algorithms are needed. Our parallel tempering cluster al- 
gorithm is a combination of replica-exchange methods f5!| 
with the Swendsen-Wang cluster algorithm 6] and is 
therefore predestinated for use in such devices. 



II. PARALLEL-TEMPERING CLUSTER 
ALGORITHM 

For the parallel-tempering procedure of the combined 
algorithm we use a set of iVrep replica, where the number 
of replica depends on the "desired" range that is needed 
for the FSS analysis . To determine this range we per- 
form at the beginning of our simulations a short run in a 



reasonable temperature interval. We choose the number 
of replica TVrep such that the acceptance rate ^(1 — >■ 2) 
between adjacent replica is about 50%, which can be cal- 
culated from 

A{1 ^ 2) = ^ Pp,{Ei)Pp,{E2)PMEi,Pi ^ i?2,/32), 

Ei,E-2 

. . . (1) 

where Pp^{Ei) is the probability for replica i at inverse 
temperature Pi to have energy Ei and 

PMEi.lii -^ E2,(32) = min[l,exp(A/3A^)] (2) 

with A/3 = 132 - h, AE = E2 - Ei is the proba- 
bility to accept a proposed exchange of different, usu- 
ally adjacent replica. This choice of A{1 — >■ 2) ensures 
that the multi-histogram reweighting [8| works properly 
and the fiow in (inverse) temperature space, that is, the 
rate of round trips between low and high temperatures 
is optimal [7|. This is the main difference to our pre- 
liminary note [9j, where we used histogram overlaps in- 
stead. Using the data of this short run as input for the 
multi-histogram reweighting routine, we determine the 
pseudo-critical points C™'^^ = C{13™'^'^) of the specific 



heat C(^) = /3^V{{e^) - 
tibility x{(3) = l3V{{m^) 
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(e)^) and ^"^^^ of the suscep- 
~ (|m|)^), where e = E/V is 
the energy density, m — M/V the magnetization den- 
sity, and V = L''' the size of the system. Furthermore, 
we measured the maxima of the slopes of the magnetic 
cumulants, C/2fc(/3) = 1 - (m2'=)/3(|m|'=)2 (fc = 1,2), 
and of the derivatives of the magnetization, d{\m\) / d(3 , 
d{\n\m\'')/d/3 {k = 1,2), respectively. We also include 
the first structure factor Sk^ (see, e.g., Ref. [10]) in our 
measurement scheme to allow a direct comparison with 
the results of Ref. [2] ■ 

Next we determine (3 values where the observables 
S" = {C, X, ■ ■ ■ } reach the value S{l3p~) = rS'^'''' with 
r < 1. This leads to a sequence of /^g values satis- 
fying /3s < /Sg"''^, and /3^ > /Sg"^'^ as illustrated for the 
two-dimensional (2D) Ising model in Fig. [T] The actual 
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FIG. 1: (Color online) Reweighted observables for the 2D 
Ising model with L = 8. The symbols mark the maximum 
values 5"""-'' and the value S{l3p~) = rS'"'''' with r = 2/3. 



simulation range is then given by the largest interval cov- 
ered by these f3g values, i.e., for the example in Fig. [T] 
the "desired" simulation window would be [/3g ,/3j]- 

In practice we start by estimating for a very small sys- 
tem a reasonable (inverse) temperature interval [/3^, /3+] 
and the number of replica A^rcp by trial and error. For 
successively larger systems we use the measured temper- 
ature interval and N^cp of the next smaller system as 
input parameters. The work flow of our method is then 
given by the following general recipe: 



1. compute the simulation temperatures of the replica 
equidistant in /3, 



2. perform several hundred thermalization sweeps and 
a short measurement run. 



3. check the histogram overlap between adjacent 
replica: if the overlap is too small (< 10%), add 
a further replica and go back to step [TJ else go on, 

4. use multi-histogram reweighting to determine 
f3g and /3^ for all observables S, leading 
to the temperature interval [/^min'Anax] — 



[m.ms{(3g},ma,^s{l3j}], 



5. start with /S = /?,-f^j,J and compute a sequence of 
temperatures /3i with fixed acceptance rate ^(1 — > 
2) until ft = /?+ > /3+,„ 

6. perform several thousand thermalization sweeps 
and a long measurement run. 



III. RESULTS 
A. 2D Ising model 

Applying this recipe to the 2D Ising model, our com- 
puter program simulated system sizes from L = 8 up to 
L — 1024 fully automatically. This shows how robust our 
new method is. Table U gives an overview of the resulting 
temperature intervals and the numbers of replica needed 
when for comparison with Ref. [31 r = 2/3 is used. Due to 
the acceptance rate criterion in step 5 of our iterative pro- 
cedure, the upper bounds /3+ of the temperature intervals 
slightly overshoot /^^ax ^^'^ show relatively large fluctu- 
ations. With increasing system sizes this discretization 
effect becomes less pronounced, see Fig. [5] where we com- 
pare the automatically determined interval of our algo- 
rithm with the exact temperature interval [f3^, f3'^] using 
the specific- heat formula of Ferdinand and Fisher [ll| . 

To assess the performance of the method, we mea- 
sured the integrated autocorrelation times Tint (ft, L) and 
determined the maximum over all replica, Tint{L) = 
maxi=ijv^^ Tint (ft, L), for each lattice size. By fitting the 
critical slowing down ansatz Tint (L) ex L^ to the data, we 
find for e, rn^, and Sk^ rather small dynamical critical 
exponents z = 0.19(1), z = 0.11(1), and z = 0.01(1), re- 
spectively. As an example, we compare in Fig.[3]our data 
for Tint,E{L) of the energy with the results of Ref. 2). Of 
course, in our case the computational effort depends lin- 
early on the number of replica needed. We therefore also 



show the effective autocorrelation time Tcff 



Nr, 



I "Hnt : 



which enables a fair comparison in units of lattice sweeps. 
We see that for L > 100, our Tcff is more than one order 
of magnitude smaller than using the recently proposed 



multibondic Wang-Landau method [4I. For r = 2/3, N^ 



rep 



grows with increasing lattice size ex L^ with z' w 0.18 
(cf. Table[l|. Consequently, Toff ex L''"" with Zoff = z + z' . 
For the energy this gives z^g « 0.37, which is still much 
smaller than the exponent z « 1.04 found in Ref. [2], 
so that the gain in efficiency becomes more and more 



TABLE I: Simulation windows and numbers of replica (see 
text) for the 2D Ising model simulations with r — 2/3 on L^ 
lattices. 



L 


r = Ps,^ /3+.. = (3+ 


/?+ 


iVrcp 


K% 


8 


0.191636 0.489 877 


0.565 730 


7 


6 


16 


0.320 580 0.470 826 


0.471618 
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7 


32 


0.380 475 0.460 499 


0.471 780 


9 
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64 


0.410 832 0.452 619 


0.457 003 


10 
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128 


0.425 789 0.448 345 


0.449 188 


11 
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256 


0.433 267 0.446 069 


0.446 392 


13 
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512 


0.437 042 0.443171 


0.443 674 


14 


10 


1024 


0.438 854 0.442 400 


0.442 464 


16 


10 


00 


Pc = ln(l + V2)/2 = 0.440 686 7935 . . . 
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FIG. 2: (Color online) The "desired" temperature interval for 
r = 2/3 as a function of the system size. The horizontal hue 
indicates the critical inverse temperature Pc, the other lines 
show exact results derived from Ref . [IJl . The circles indicate 
the simulation windows determined fuUy automatically, cf. 
Table |IJ and the boxes show for completeness the measured 
values for /3+, /3J?"" and P^- 
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FIG. 3: (Color online) Autocorrelation times rint and Tcfi for 
the energy of the 2D Ising model, where TcH = A^rop Tint with 
A'^i-cp denoting the number of replica, cf . Table ID 



pronounced with increasing lattice size. 

The choice r = 2/3 of Ref. [21 is quite conservative as 
even for r = 1 the peaks determining the left and right 
boundary of the "desired" simulation window are usually 
sufficiently well sampled. This amounts to a somewhat 
smaller temperature range and repeating the above pro- 
cedure with r = 1, we arrive at a smaller ZcS — z = 0.27, 
see Fig. |3l Here ZcS = z because the integer valued 
iVrop turn out be so small (3 for L — 8, 16 and 4 for 
32 < L < 1024) that they practically stay constant over 
a wide range of system sizes. For large L, this leads to a 
significant gain, and also for the moderate system sizes 
of Fig. 131 Toff is already reduced by a factor of about 3 
compared with the case of r = 2/3 and a factor of about 



FIG. 4: (Color online) FSS of the "desired" simulation win- 
dow for the 2D Ising model with r = 2/3. 



100 compared with Ref. [21 ■ 

On the other hand, for r < 2/3, the "desired" simula- 
tion window (and hence also iVrop) becomes larger, and 
the effective autocorrelation time Tcff grows faster with 
system size. For example for r = 1/2, we find z' ~ 0.32. 
With z = 0.23(1) this leads to z^s = 0.55(2) for the 
energy and similar results for m'^ and Sk-i ■ 

In order to understand the FSS of Tcb theoretically, we 
have to recall that the 2D Ising model is a particular case 
due to its logarithmic divergence of the specific heat C 
(with critical exponent a = 0). First, this leads to a log- 
arithmic correction in th e FSS of the reweighting range, 
A/^rew (X L~^l" I\J\ylL. Second, when C is included 
in the "desired" quantities, it (empirically) determines 



/3i, which then does not scale 



121. For the "de- 

i/'^+foL-''/'', this 



the upper bound (i'^^ 

generically with L~^l'" but with L~'^l" 

sired" range /3+—/3- w /3j—/3g — aL 

shows that L~'^l" is asymptotically the leading ter m for 
r < 1, so that iV^op = (/3+-/3-)/A/3rcw ^ L^i-'-)/'^^!^ 
for L > Lx = [ajhyl'^^^'^^ . However, as Fig. |3| demon- 
strates for r = 2/3, this is not observed for practically 
accessible lattice sizes: Instead of L~'^l'^ (with z^ = 1) we 
observe /S"*" — /3^ ex L^" with k' — 0.94. The reason is 
the very slow crossover from subleading to leading scal- 
ing behavior. In fact, whereas both Pc — Pg ex. L~^-°^ 

and /3^ — /S^ oc L~'^-^^ do scale as expected, the asymp- 
totically leading term starts to dominate only at around 
L = 1000 ~ Lx with a/b — 10, which implies A^rcp oc 
L^'^^VuiL for L < 1000. Since in this range to a very 
good approximation vhiL ^ L^-^^ (which can be easily 
verified directly), this effectively leads to A^rep oc L°-^'^, in 
good agreement with our direct estimate z' « 0.18. For 
r — 1/2, L~^ and L~^ = L"^/^ differ more strongly and 
we observe the crossover in the FSS of /3+ — /3^ already 
at about L = 150 w L^ with a/b ~ 12, yielding here 
P+ -p- (X L-0-«o and hence TV^op oc L^-^^Vh^L « iO-^i, 
again as fitted directly. Only for r = 1, /3+ — /?^ ex L^^/'^ 



and iVic 



'mL, which as explained above is hardly 



visible for small (integer) values of iVrop- 

If one omits C as a criterion to specify the "desired" 
FSS window, both the upper bound /3+ and /3+ — /3~ 
scale oc L^^^'^ for any value of r, so that N^ 



rep 

see the last column in Table |T] for the case r — 2/3. If 
the simulation window is now too narrow to determine 
the critical exponent a directly, one still can use the hy- 
perscaling relation a = 2 — diy with d the dimensionality. 
Whereas the "bare" dynamical critical exponents z do 
not differ much, due to the smaller number of replica 
needed, TeS = N*^^ Tint is slightly smaller than in the 
case including C. If one again simply fits a power law to 
both Tint and Tcff (i.e., ignores the logarithmic behavior 
of A^*ep), one finds excellent fits with z = 0.20(1) and 
z^ff = 0.31(1) = z + 0.11 for the energy (z = 0.11(1), 
Zeff = 0.22(1) for m^ and z = 0.03(1), z^s = 0.13(1) 
for Ski), confirming that effectively VhiL sa L^-^^ for 
moderately large L < 1000. 



B. 3D Ising model 

In the 3D Ising model where a « 0.11 > 0, both the 
reweighting range and the "desired" temperature window 
should scale with L"^/", so that one would expect that 
our routine will use the same number of replica for all 
system sizes. This is indeed the case for the simplest 
choice r = 1, where our routine determines for all lattice 
sizes A'icp — 4. The autocorrelation analysis along the 
same lines as in 2D gives z = z^s = 0.62(2) for the energy 



TABLE II; Simulation windows and numbers of replica for 
the 3D Ising model simulations with r = 2/3 on L^ lattices. 
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/5--/^s„ 


Pi 


/3+ 


-^Vrcp 
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0.061 955 


0.284 066 


0.311321 
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0.142 959 


0.259 876 


0.272 522 
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8 


0.173 712 


0.252 866 


0.259 229 
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10 


0.188 447 


0.246 916 


0.252 136 


10 


12 


0.197 404 


0.241 159 


0.243 106 


10 


16 


0.206 422 


0.236 610 


0.236 773 


11 


20 


0.211183 


0.233 407 


0.233 621 


12 
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0.213 807 


0.232 462 


0.233 016 


14 


28 


0.215 553 


0.229 853 


0.230 301 


14 


32 


0.216 670 


0.229 221 


0.229 293 


15 


36 


0.217 572 


0.228 026 


0.228 613 


16 


40 


0.218 172 


0.227 703 


0.228 025 


17 


48 


0.219 082 


0.226 669 


0.226 754 


18 


56 


0.219 621 


0.225 353 


0.225 555 


18 


64 


0.220 031 


0.224 758 


0.224 775 


18 
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0.220 309 


0.224 359 


0.224 435 


19 


80 


0.220 505 


0.224 331 


0.224 347 
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13c ~ 0.221 654 59 (Ref. fl3]) 
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FIG. 5: (Color online) Same as Fig.|3]for the 3D Ising model, 
cf. Table nil 



{z = zoff = 0.59(2) for m^ and z = z^s ^ 0.37(2) for S'fcJ. 
This exponent is thus again much smaller than z « 1.05 
obtained in Ref. [2|, and already for moderate system 
sizes L « 40 — 80 the values of Tcb are about 20 — 30 
times smaller, cf. Fig. [5] 

If we choose r = 2/3 as in Ref. [2|, however, we find also 
here a weak system-size dependence A'rep oc L^-^^, cf. Ta- 
ble [Tl] where also the automatically determined tempera- 
ture intervals are given. Here we obtain z = 0.44(1) and 
thus Zeff = 0.80(1) for the energy, cf. Fig.[5](z = 0.41(1), 
z^g = 0.78(1) for m2 and z = 0.18(2), z^s = 0.55(2) 
for Ski ) • The reason for this unexpected result can be 
traced back to the fact that the upper boundary /3p of 
the "desired" simulation window, determined by the low- 
temperature tail of C, lies for r ~ 2/3 clearly outside the 
FSS region. In fact, omitting C as a criterion for the FSS 
window, iVj-cp = 7 — 9 stays almost constant. The choice 
r < 1 is thus less favorable, but even for r — 2/3 one 
gains about one order of magnitude in computing time 
compared with Ref. [2,]. 



IV. CONCLUSIONS 

To summarize, we have introduced a very flexible and 
simple approach for a systematic determination and sim- 
ulation of the critical temperature window of interest 
for second-order phase transitions, which one needs for 
an accurate estimation of critical exponents and other 
quantities characterizing critical phenomena. The effi- 
ciency of the method depends of course on the chosen 
or available update scheme in the particular case, with 
non-local cluster flips being the favorable choice. Since 
the setup of our method is completely general and can be 
combined also with any other update scheme (multigrid, 
worms. Metropolis, heat-bath, Glauber, . . . ), it could be 
employed for all simulations in statistical physics, chem- 
istry and biology, high-energy physics and quantum field 
theory where one is interested in critical phenomena. 
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